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Abstract 

We have designed a new symbolic-numeric strategy to compute efficiently and accurately floating 
point Puiseux series defined by a bivariate polynomial over an algebraic number field. In essence, 
computations modulo a well chosen prime p are used to obtain the exact information required to 
guide floating point computations. In this paper, we detail the symbolic part of our algorithm : 
First of all, we study modular reduction of Puiseux series and give a good reduction criterion to 
ensure that the information required by the numerical part is preserved. To establish our results, 
we introduce a simple modification of classical Newton polygons, that we call "generic Newton 
polygons", which happen to be very convenient. Then, we estimate the arithmetic complexity 
of computing Puiseux series over finite fields and improve known bounds. Finally, we give bit- 
complexity bounds for deterministic and randomized versions of the symbolic part. The details 
of the numerical part will be described in a forthcoming paper. The reader is referred for now 
to a preliminary version sketched in (Poteaux (2007)). 
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1. Introduction 

Let K be a number field and F(X, Y) be a bivariate polynomial in K[X, Y] such that : 

• F(X, Y) = Y d + A d - 1 (X)Y d - 1 + ■■■ + A (X), with d > and n = deg x F > 0. We 
denote by D the total degree of F. 

• F is squarefree. 

• Ap(X) denotes the discriminant of F with respect to Y . 
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A root of Ap will be called a critical point. Critical points can also be defined as the set 
of numbers xo such that F(xq, Y) is not squarefree. The equation F(X, Y) = defines d 
algebraic functions of the variable X , which are analytic in any simply connected domain 
T> C C free of critical points. 

Values of these algebraic functions can be computed by means of analytic continuation 

(see Chudnovskv and Chudnoyskvl |l986l ). ICnudnovskv and Chudnovskv ( 1987 ) lvan der Hoeven 
(1999), Ivan der Hoevenl (|2005l ). who require that a differential equation satisfied by the 
functions is known). In our context (see below), no differential equations is known a priori 
and we do not need a high precision since we just need enough information to separate 
functions. Hence, it is not clear that these asymptotically fast (in terms of precision) 
methods are relevant. 

However, if T> is included in a sufficiently small disc centered at a critical point xo, it 
is well-known that numerical values of these functions in T> can be obtained directly via 
truncated Puiseux series at X — xo (see Section 2). 

We have used this fact to devise an algorithm to com pute the monod romy of the Rie- 
mann sphere covering defined by the curve F(X, Y) = (iPoteauxl (|2007l) V The algorithm 
follows paths along a minimal spanning tree for the set of critical points and expansions 
above critical point are used to b ypass them. 

In (|Deconinck and van Hoeii (2001)), a different strategy for computing the mon- 
odromy is described. The approach therein and the implementation suffer from numer- 
ical instability and a lack of error control. These weaknesses have stimulated us for 
investigating the method mentioned above, that could hopefully offer a better numer- 
ical control. Another application of Puiseux series is the computati on of Abel's m ap, 
at the heart o f the celebrated Abel- Jacobi theorem (see for instance iMirandal (|1995l ) or 
Forster ( 198ll )). Abel's map yields particularly interesting solu tion s of certain differential 



equat ions from Physics (see iDeconinck and Patterson ( 2007 ) and Deconinck and Segur 
(1993)), and allows to answer questions about linear equivalence of divisors on algebraic 
curves. The latter has application to C ompute r Alge b ra, notably to th e determination 
of antiderivatives of algebraic functions (|Traeerl (|l984l ). lBronsteinl dl990h'). or the charac- 
teriza ti on of algebraic solutions of l inear differential equations (jBaldassarri and Dwork 
( 1979f ). Compoint and Singerl ( 19981 )1. In the context of Abel's map, one is led to com- 
pute integrals of algebraic functions along paths ending at cri tical points. Series expan- 
sions a bove critical points are definitely useful for this task (see IDeconinck and Patterson 
(2007)). Moreover, representing algebraic functions as truncated Puiseux scries could lead 
to error bounds for the integrals. 

Unfortunately, applying a floating point Newton-Puiseux algorithm (see Section 3) to 
compute Puiseux series above a critical point is doomed to failure. Indeed, if the critical 
point xq is replaced with an approximation, expansion algorithms return approximate 
series with very small convergence discs and do not retain important information, such 
as ramification indices. Therefore, the output is not aidful. 

On the other hand, coefficient growth considerably slows down symbolic methods. 
Since the degree of Ap is in 0{D 2 ), a critical point xq may be an algebraic number with 
large degree. Puiseux series coefficients above xo belong to a finite extension of K whose 
degree over K may be in 0(D 3 ). For D = 10, the extension may already be excessively 
large. Moreover, when these coefficients are expressed as linear combinations over Q, the 
size of the rational numbers involved may also be overwhelming. Floating point evaluation 
of such coefficients must, in some cases, be performed with a high number of digits 
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because spectacular numerical cancellations occurs (see examples in IPoteauxl ( 2007l )l. 
For instance, the polynomial of degree 6 in Y F(X, Y) = (Y 3 - X) ({Y - l) 2 — X) (Y — 
2 - X 2 ) + X 2 y 5 has a discriminant A F (X) = X 3 P(X), where P(X) is an irreducible 
polynomial of degree 23 over Q. Rational Puiseux series (see Section 2) above P(X) have 
cocfhcients in a degree 23 extension of Q. Rational numbers with 136 digits appear in 
the first term of the expansions. Walsh has shown that, for any e > 0, the singul ar part 
of Puiseux series can be computed using 0(d' i2+f n 4+c log /i 2+c ) bit operations (Walsh 
4200(A ). where h is the height of F. Although this bound is probably not sharp, it is not 
encouraging and tends to confirm experimental observations. 

To alleviate these problems, we introduce a symbolic-numeric approach : exact impor- 
tant information is first obtained by means of computation modulo a well chosen prime 
number p, then this information is used to guide floating point computation. The coef- 
ficient size is therefore kept under control while numerical instability is reduced. Exact 
important data, such as ramification indices and i ntersection m ultiplicities of branches, 
are preserved. Experimental evidences reported in ([Poteauxl (|2007l )l seem to validate this 
approach. 

This paper presents several contributions : 

• We introduce "generic Newton polygons" and "polygon trees" (Section 3). The latter 
capture precisely the symbolic information needed for floating point computations. 

• We study modular reduction of Puiseux series and rational Puiseux expansions. This 
leads to a fully proved and easy to check criteria for the choice of a "good prime" p 
such that polygon trees can be obtained using modular arithmetic (Section 4). We rely 
on technical results that are proven or recalled in Section 2 and 3. 

• Complexity estimates for the computation of singular parts of Puiseux series over finite 
fields are given in Section 5. They improve previously published results. 

• Finally, we study the bit-complexity of a Monte Carlo approach for the symbolic part 
of our symbolic-numeric method in Section 6. 

Obtaining floating point Puiseux series from po lygon trees is n ot a trivial task. The 
principles of our technique are briefly described in ( Poteaux ( 20071 )1. Several variants, as 
well as round-off error control and overall complexity, are still under investigation. These 
questions are left to a forthc oming paper. Moreover, an implementation improving the 
prototype of ( Poteaux ( 20071 )') is being developed. 

Modular methods are extensively used in Computer Algebra to avoid intermediate 
coefficient swell (via Hensel li f ting o r the Chinese Remainder Theorem, see for instance 
von zur Gathen and Gerhardl l| 19991) 1. However, a "reduce mod p and lift" or "reduce 



mod pi and combine" method would not help much in this case since we are not facing 
an intermediate coefficient growth problem, but an intrinsically large symbolic output. 
Modular methods are much less common when it comes to directly obtaining numerical 
results. We are aware of an attempt by Jean-Charles Faugere to control Grobner bases 
computations that does not seem to have been successful so far (personal communica- 
tion). We therefore claim some originality with this approach. 

To conclude this introduction, we introduce notations and assumptions that will be 
used throughout the paper. We also recall well-known facts : 

• All fields considered in the paper are commutative. 

• If L is a field, L will denote an algebraic closure of L. 

• For each positive integer e, Ce is a primitive e-th root of unity in L. Primitive roots 
are chosen so that C^ b = ( a . 
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• vx denotes the X-adic valuation of the fractional power series field L((X 1 / 6 )), nor- 
malized with vx{X) = 1. If S e L^X 1 ^)), we denote by tc(5) the trailing coefficient 
of S 1 , namely S 1 = tc(S)X Vx< - s ^ + terms of higher order. 

• If S = J\ ®-kX k l e is a fractional power series in L((X 1 / 6 )) and r is a rational number, 
S r denotes the truncated power series S r = J2k a kX k/e where TV = max{/c e N | f < 
r}. We generalize this notation to elements of L^X 1 ^))^} by applying it coefficient- 
wise. In particular, if H e L[LY]][Y] is defined as H = J2i(J2k>o a ikX k )Y l , then 

• The discriminant of a univariate polynomial U is denoted by Ay. If U is a multivariate 
polynomial, the context will always allow to identify the variable. 

• If U(T) (resp. V(T)) is a monic separable univariate polynomials of degree s (resp. r) 
with roots {u\, . . . , u s } (resp. {vi . . . , v r }), then : 

A n = [~| (m-Uj) Resultant (U, V) = ]] (m-Vj). (1) 

1<»,J<S l<i<s 
i^j l<j<r 

• If U is a monic univariate polynomial which admits a factorization into a product of 
monic polynomials U — 111= 1 ^ then : 

r 

Al/ = ri A ^ II Resultant (Ui, C/,). (2) 

i=l l<i,J<r 

• For each i with < i < e — 1, we denote [i, e] the automorphism : 

M : L((XV«)) ^ I((xVe)) 

When the ramification index can be deduced from the context, we shall simply write 
[i] instead of [i,e]. If S 1 e L{{X^ e ), the image of S under [i] is_denoted by S®. This 
notation extends naturally to any polynomial with coefficient in L((X 1 / e ). It is obvious 
that elements of the subficld L((X)) are invariant under [i]. 

• Let / be a polynomial in L[T] with squarefree factorization / = Yii=i fi l (that is, the 
ki are pairwise distinct positive integers and the fa are positive degree polynomials 
with gcd(/j,/j) = 1 if i ^ j). We associate to / the partition of deg/ denoted 
[/] = (/cf° s ^ . . . fcj? cg ir ) . Namely, the multiplicity ki is repeated deg fi times in the 
decomposition of deg/. 

• If H £ L[X, Y], then Hx and Hy are the formal partial derivatives of H. 

• For a multivariate polynomial H(X_) = J2 k a k 2L~ € C\X] = C[X±, . . . , X n ], where k is 
a multi-index, we denote : 



\H\\ 00 =m*x{\a k \} \\H\\ 2 = JY,{\^\ 2 - 



2. Puiseux series 

We need to state results over more general fields than K. Throughout the section, L 
stands for a field of characteristic p > and F belongs to L[X, Y]. Otherwise, we keep 
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the assumptions and notations of Section 1. We also impose the condition : 

p = or p> d = dcg Y (F) (3) 
Up to a change of variable X <— X + x n , we assume that the critical point is X = 0. 

2.1. Classical Puiseux series 

In this part, we review classical results about Puiseux series. We begin with : 

Theorem 1 (Puiseux). Let H be squarefree polynomial of L[X,Y] such that deg Y (H) = 
d>0. 

• If condition (3) is satisfied, there exist positive integers ei, . . . , e s satisfying X)i=i e i = d 
such that H (viewed as a polynomial inY) has d distinct roots in L((X)) which can 
be written : 

s ij (x)= ^(t-x^ 

k> — oo 

for 1 < i < s and < j < ej — 1. Moreover, the set of coefficients {aik} is included in 
a finite algebraic extension of L. 

• If p = 0, then : 

ifflclp))= (J I((* 1/e )) 

eGN* 

Definition 2. These d fractional Laurent series are called Puiseux series of H above 
0. The integer is the ramification index of Sij. If > 1, then 5V, is ramified. If 
Sij G L^X 1 /^]], we say that SVj is defined at X = 0. If SVj(O) = 0, we say that Sij 
vanishes at X = 0. 

An arbitrary number of terms of all Puiseux series can be effectively computed using 
Newton-Puiseux algorithm (see Section 3). 

For each positive integer e < d, our hypothesis (3) about the characteristic p of L imply 
that the Galois group G e of L((X 1 /^))/L((X)) is cyclic and generated by [1] : X x l e ^ 
CeX 1 ^. Hence, G ei permutes cyclically the elements of the set 5^ = {Sij(X)} <j< ei -i. 

Definition 3. We call Si a cycle of H above 0. If an element of Si (thus, all elements) 
vanishes at X = 0, we say that the cycle vanishes at X = 0. 

Since the SVj (0 < j < — 1) can be quickly recovered (both symbolically and 
numerically) from any element of Si, it is sufficient for our purposes to compute a set of 
representatives for the cycles of H. 

Definition 4. The regularity index of Sij in H is the least integer N such that 

. JV JV 

Sij" i — S uv " i implies (u, v) — (i,j). The truncated series SYj ei is called the singular 
part of Sij in H. 

In other words, rij is the smallest number of terms necessary to distinguish Sij from 
the other Puiseux series above 0. It is worth noting that depends not only on Sij, but 
also on H since H is not assumed irreducible in L[X, Y\. See examples in Section 3.2. 
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If the singular part of a Puiseux series is known, a change of variable yields a bi- 
variate polynomial for which re maining terms of the ser ie s can be computed "fast" us- 
i ng qu adratic Newton iterations (|Kung and Traubl (|l978l ). Ivon zur Gathen and Gerhardl 



(|1999T )). Newton iterations can be applied to series with floating point coefficients, there- 
fore we focus on the computation of the singular parts of the Sij . Since it can be shown 
that all elements of a cycle Si have the same regularity index, that we denote r,, the 
problem reduces to the determination of the singular part of a representative of Si for 
1 < i < s. 

Since F is monic, all Sij arc defined at X = 0. In other words, they are fractional 
power series. In this case, when L C C, the Sij have a convergence radius which is at 
least equal to the distance from to the nearest (nonzero) critical point. If we choose 
a determination for the e^-th root functions, the Sij define d analytic functions in any 
domain T> that is included in the convergence disc and does not intersect the branch cut. 
To evaluate accurately these functions in T> 7 we need to : 

• Control truncation orders of Puiseux series. Bounds are given in (|Poteauxl (|2007h ). 



• Compute efficiently floating point approximation of th e truncat e d Sg . The principles 
and a preliminary implementation are presented in (|Poteauxl (|2007t n. The present 
paper details the symbolic part of our method. 

• Give error bounds for the approximations of the aik and study the algorithm numerical 
stability. This topic has not been addressed yet. 

2.2. The characteristic of a Puiseux series 

We derive relations between the discriminant of F and particular coefficients of its 
Puiseux scries that we shall use to define a "good reduction" criterion. Let 

oo 

s(x) = J2^ l/e 

denote a Puiseux series of F with ramification index e > 1. We define a finite sequence 
(B ,R ), (Bi, . . . , (Bg,R g )) of integer pairs as follows : 

• R = e, B = 0. 

• If Rj-i > 1, we set Bj = min {i > Bj-\ | a, 7^ andi ^ (mod Rj — gcd(_Bj, Rj-i) 
If Rj-i = 1, we stop and set g = j — 1. Note that g > 1 and R g = 1. 

Finally, we set Qj = Rj^i/Rj, Mj = Bj/Rj for (1 < j < g) and define Hj to be the 
largest non negative integer so that Bj + HjRj < Bj + i for < j < g — 1. It is clear that 
e = Q1Q2 ■ ■ ■ Q g and Mj is an integer prime to Qj. 

With these notations, and up to a new indexing of the coefficients, S can be written 
in the form : 

S(X) = E§o«o^ 

"l „ 77, M l+1 

M 2 it M 2 +j 



+ l2 x^ + a 2d X^ (4) 



Mg +j 



In the expression above, the monomials of S are ordered by increasing (rational) 
degree. 



G 



Definition 5 ((see Zariskil (|l981 ) or Brieskorn and Kndrrer ( 1986h ). The characteristic 
of S is the tuple of integers (e; f?i, . . . , B g ). The characteristic coefficients are the ele- 
ments of the sequence (71, . . . ,j g ) and characteristic monomials are the corresponding 
monomials of S. 



Proposition 6. Assume that hypothesis (3) is satisfied. Let G(X,Y) be the minimal 
polynomial over L((X)) of a ramified Puiseux series S G i[[X 1/e ]] as above. Let A G (X) 
be the discriminant of G with respect to Y . Then : 



tc(A G )=±(nQf i n^ i " fl< ) e 



i=l 1=1 



(5) 
(G) 



Proof. We first introduce the notation v = vx(Aq) ancl — tc(A G ). The conjugates of 
S over L((X)) are {5 , ^}o<i< e -i, therefore : 

A G = [] (S®-SW). 

0<i,j<e-l 

From this relation and (5), we note that v depends only on the contribution of the 
terms X Bi ^ e — X Mi < W 1 '" 1 ^), Hence, v is determined by the exponent of ji in 9. There- 
fore, if (5) is true, so is (6) since : 

a B 9 

v = J2 e CRi-i - R >)^T = Yl B i( R i-i ~ Ri)- 

i=l i=l 

In order to prove (5), we proceed by induction on g. For each positive integer r let S r 
be the discriminant of X r — 1, that is S r = ±r r . 

If g = 1 , the expansion of Aq in increasing fractional power of X is : 



a g = n (j g (c^~c^)x MjQ °+---) 



0<i,j<c-l 



s(e-l) 



n cc 

, 0<»,3<e-l 



Mr, 



x e(e-l)M g /Q g _, 



Since M g is prime to Q g and Q g — e, is a primitive e-th root of unity. We obtain 
9 = 8 el e g (e - 1] = ±Qf 1 e g {R *- 1 ~ R * ) as expected. 

We now assume that g > 1. To simplify notations, we set Q = Qi and R = Ri = 
Q 2 • • • Q 9 - We define if G LffX 1 /* 3 ]]^] as follows : 

H-l 



if = J](y-5 [iQ] ). 



i=0 
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Since [Q] = [Q,e] generates the Galois group of L^X 1 ^)) over L^X 1 ^)), H is the min- 
imal polynomial of S over L((X 1 / < 3)). Moreover, the factorization of G over L((X 1/,( 2)) 
is given by : 

Q-i 

g= n ^w. 

z=0 

Using relation (2), we obtain Aq = nili2 where : 
Q-i 

Hi=JjA HW II 2 = JJ Resultant ii^l). 

i=0 o<i,j<Q-i 

We need to evaluate the contribution to 9 of III and n 2 . We first consider EEi. Let 
U(X,Y) = H{XQ,Y) be the minimal polynomial of S(X Q ) over L((X)). Since U has 
characteristic (i?; P>2, ■ ■ ■ , B g ), our induction hypothesis yields : 

i=2 i=2 

for some positive integer u. Therefore : 

A Hl3] =C R3 (±f[Q«' f[7« i - 1 - Ri ) R X%+... 

j=2 i=2 

Since Q R = e and C™ ^ = Cq J the contribution of III to 9 is : 

j=2 i=2 

We now estimate the contribution of Resultant (i? W, ff^). Each difference of roots in 
the product defining the resultant has the form : 

7i(Cr-C^')(^ Ml/gi + ...) 

and there are R 2 such differences. Since there are Q(Q— 1) resultants in the product and 
i^f 1 is a primitive Q-th root of unity, we conclude that the contribution of II2 to 9 is : 

D<«,3<Q-1 

Combining the last expression with (7) gives (5). □ 

Remark 7. The value of ux(Aj?) is well-known (see for instance ( Zariski ( 198li 0. It can 
be expressed as the sum of a differential exponent and of a conductor degree. In Singular- 
ity Theory, it also has an interpretation in terms of "infinitely near point" multiplicities. 
However, we have not found in the literature an expression for tc(Af ). 

2.3. Rational Puiseux expansions 

In order to perform computations in the smallest possible extension of L and to take 
advantage of c qnjugacy over L, Duval introduced the notion of "rational Puiseux expan- 
sions over L" (|Duva]| (|l987h ). This arithmetical concept is irrelevant in the context of 
floating point computations, but will prove useful for expansions over finite fields. 
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Remark 8. A diffe rent de f initio n of "rational Puiseux expansions over V appeared in 
(iDuvall (|l989h ) and (IWalsbl ( 199Sh ). The definition given therein corresponds to "rational 



Puiseux expansions over L" in the sense of (|DuvaJ (|l987l) l and in the sense of this paper 



Definition 9. Let H be a polynomial in L[X, Y] with degy H > 0. A parametrization 
R(T) of H is a pair of non constant power series R(T) = (X(T),Y(T)) E L((T)) 2 such 
that H(X(T),Y(T)) = in L((T)). The parametrization is irreducible if there is no 
integer u > 1 such that i?(T) e L((T")) 2 . The coefficient field of -R(T) is the extension 
of L generated by the coefficients of X(T) and Y(T). 

Assume for a moment that H is irreducible in L[JT,y] so that JC = L(X)[Y]/ (H) 
is an algebraic function field. A parametrization R(T) = (X(T),Y(T)) induces a field 
morphism : 

<P R : JC I((T)) 

/(X,Y) i ^ f(X(T),Y(T)) 

Composing with the valuation t>*r of L((T)), we obtain a valuation of JC that we denote 
again by Vt- It is easily seen that the set = {/ G /C | ^t(/) > 0} is a pZace of /C in 
the sense of dChevallevi (<195ll) 1 and that Vr = {/ G /C | v(f) > 0} is the corresponding 



V-ring of JC. We recall that is the unique maximal ideal of Vp and that the residue 
field of is Vr/^Pa, which can be viewed as a finite algebraic extension of L. Therefore, 
we obtain a mapping ^ from the set of parametrizations of F onto the set of places of 
JC. Reciprocally, to each place of JC corresponds parametrizations of H . 

Let us denote by {^5,:}i<i< r the places of JC dividing X and by ki the residue field of 
%■ 

Definition 10 (Rational Puiseux expansions). 

• Assume that H is irreducible in L[X, Y]. A system of L-rational Puiseux expansions 
above of H is a set of irreducible parametrizations {Ri}i<i< r of the form : 

R i {T) = {X i {T),Y i {T))= UT% Y, ^ Tk ) ^{{T)f 

with e,i > such that : 
(i) ^ is one-to-one from {Ri}i<i< r to {tyi}i<i< r . We assume that the ^ are numbered 

so that q3i = =^(Ri). 
(ii) The coefficient field of Ri is isomorphic to ki. 

• Assume that H is squarefree. A system of L-rational Puiseux expansions above of H 
is the union of systems of L-rational Puiseux expansions for the irreducible factors of 
H in L[X, Y]. 

Definition 11. We say that Ri is defined atT = \iYi £ L[\T]]. In this case, the center 
of Ri is the pair (X t (0), r,(0))el 2 . 

The classical formula rel ating degrees of r esidue fields and ramification indices of an 
algebraic function field (see Chevallev ( 195 if )) translates into : 
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Theorem 12. Let H e L[X, Y] be squarefree and deg Y H = d > 0. Let {Ri}i<i< r be a 
system of L-rational Puiseux expansion above for H . Let fi stand for the degree over 
L of the coefficient field of Ri . Then : 



i=l 

Classical Puiseux series can be readily deduced from a system of rational Puiseux 
expansions : 

(1) Ri has exactly fi conjugates over L, that we denote R? (1 < a < fi). 



(2) Each Rf yields a Puiseux series Yf((X/j" ) 1/<ei )- The set of all such series form a 
set of representatives for set of cycles {Si}i<i< s of H above 0. 

(3) The d Puiseux series are finally obtained using the action of G ei , 1 < i < s. 

In particular, classical Puiseux series that are defined at X = (resp. that vanish at 
X = 0) correspond to rational Puiseux expansions defined at T = (resp. centered at 



Again, we note that regularity indices for all Puiseux series corresponding to the same 
rational Puiseux expansion are equal. Therefore, we define the singular part of a rational 
Puiseux expansion Ri to be the pair : 



where ri is the regularity index of a Puiseux series associated to Ri . 

Finally, we note that since F is monic, rational Puiseux expansions of F above zero 
are all defined at T = 0, that is Ri(T) eL[[T}} 2 . 

3. Symbolic algorithms for Puiseux series 

In this section, we describe two algorithms to compute singular parts of Puiseux 
series and rational Puiseux expansions. Both methods are used by our symbolic-numeric 
approach. We also explain how coefficients computed by the two methods are related and 
conclude by comments about coefficient growth. 

Throughout the section, L stands again for a field of characteristic p > and F is a 
polynomial L[X, Y] such that condition (3) is satisfied. Moreover, we keep the assump- 
tions and notations of Section 1. 

Newton polygons and characteristic polynomials are the crucial tools. We first re- 
call well-known definitions and introduce a variant that will prove more convenient and 
powerful. 



r 





(0,0)). 
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3.1. Generic Newton polygons and characteristic polynomials 

Assume that H(X, Y) = J2 hj a^X^Y 1 is a polynomial of L[[X]] [Y] such that H (0, Y) ^ 
0. The Newton polygon of H is classically defined as follow : 

Definition 13. Denote by 1(H) the nonnegative integer vy(H(0,Y)) and by TL the 
convex hull of Supp(#) = {(i, j) £ N 2 | a i3 ^ 0}. The Newton Polygon Kf(H) of H is the 
lower part of TL. Namely : 

• If H(X, 0) 7^ 0, N(H) is formed by the sequence of edges of TL closest to the origin 
and joining (0, v x (H(X, 0))) to (1(H), 0), 

• If 0) = 0, (0, vx(H (X, 0))) is replaced by the leftmost point of TL with smallest 
j-coordinate. 

In particular, a vertical or horizontal edge is not considered part ofAf(H). The Newton 
polygon may be reduced to a single point. For instance H(X,Y) = Y yields the trivial 
polygon (1,0). 

We now introduce a slightly different object, that we call generic Newton polygon for 
reasons explained later. This variation allows a homogeneous treatment of finite series, 
clearer specifications for the algorithms and simplifies wording and proofs of results 
regarding modular reduction. 

Definition 14. The generic Newton polygon QN(H) is obtained by restricting N(H) 
to edges with slope no less than —1 and by joining the leftmost remaining point to the 
vertical axis with an edge of slope —1. 

In other words, we add a fictitious point (0, jo) to Supp(_ff) so as to mask edges with 
slope less than -1. 

Example 15. Consider Hi(X, Y) = Y 7 +X 2 Y 2 +XY 4 +X 8 +X 6 +Y 3 X 2 +Y 5 X+Y 3 X 4 . 
In Figure 3.1, the support of Hi is represented by crosses, QAf(H\) is drawn with plain 
lines while the masked edge of M(Hi) is represented by a dotted line. 

/is 

X 



6X 




Fig. 1. £./V(#i) versus N(H{) 
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Example 16. Consider H 2 (X,Y) = Y 8 + 3X 2 Y 3 + XY 5 + YX 8 + 2X 6 + AY 2 X 3 + 
Y 5 X 2 + Y 5 X 3 + Y 3 X 4 and Figure 3.1. The edge with slope -1 is prolongated until the 
vertical axis. 

X 



X 




Fig. 2. QN{H 2 ) versus N{H 2 ) 



Example 17. Assume that Ha(X,Y) — Y. Then QN(Ha) is formed of a unique edge 
joining (0, 1) to (1,0). 

The first stage of the algorithms requires a special treatment. To this effect, it is 
convenient to introduce the following definition : 

Definition 18. The exceptional Newton polygon £M(H) is the unique horizontal edge 
[(0,0),(deg y (ff(0,r)),0)]. 

In particular, £Af(F) = [(0, 0), (d, 0)] since F is monic. 

To an edge A of QAf(H) (or Af(H), £J\f{H)) corresponds three nonnegative integers 
q, m and I with q and m coprime such that A is on the line q j + m i — I. If A is the 
horizontal edge of £N(H) , m = I = and we choose q = 1 . 

Definition 19. We define the characteristic polynomial 0a : 

<M(T) = ^2 a ij T ^ 
(i,j)eA 

where io is the smallest value of i such that belongs to A. 

Note that if M{H) is used, (/)a{T) cannot vanish at T = 0, while QJV(H) allows 
such cancellation if A is a fictitious edge (or contain a fictitious part). In this case, the 
multiplicity of as a root of <Pa (T) is the length of the fictitious edge (or portion of edge) 
added. For £J\f(H), can also be a root of the characteristic polynomial. 

The next two lemmas recall the relation between Newton polygons of H and Newton 
polygons of its factors in L[[X]][Y] : 

Lemma 20. If H is an irreducible polynomial of L[[X]]\Y] and H(0, 0) = 0, then QM{H) 
has a unique edge A and 4>a has a unique root. 
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Proof. This is well-known for classical Newton polygons (jBrieskorn and Knorrei 
The extension to generic polygons is straightforward. □ 



(1986)). 



Lemma 21. Let Hi and H2 be elements of L[[X]][Y]. Then, QN^HxH-z) results from 
joining together the different edges ofQM{H\) andQN(H2), suitably displaced. Moreover, 
the characteristic polynomial of an edge A with slope —m/q of QAf(H±H2) is the product 
of the characteristic polynomials associated with edges of QN{Hi) and QN{H2) of slope 
—m/q. In particular, if Hi(0,0) ^ 0, so that GN{H\) is reduced to the point (0,0), then 
gN(H 1 H 2 )=QN(H 2 ). 



Proof. For classical Newton polygons, see ([Brieskorn and Kndrrer (1986)). For generic 



Newton polygons, proceed as follow : If necessary, add a monomial cX Ul (resp. cX U2 ) to 
Hi, where c is an indeterminate, so that QN{Hi) = jV(flj). Then, apply the result for 
the classical case and set c = to recover QAf(HiH2). □ 

Both algorithms below perform successive changes of variable, determined by (q, m, I) 
and the roots of 0a- They return a set of triplets {(Gi(X,Y), P i (X) 1 Q i (X,Y))} i such 
that : 

• Gi,Pi,Qi e L[X,Y], 

• Pi{X) is a monomial of the form AjX ei , 

• Qi(X,Y) — Qio(X) + Y X Vi , where is the regularity index of the expansion and 
(Pj(T), Qio{T)) is the singular part of a parametrization of F, 

• There exist nonnegative integers Li such that Gi(X, Y) = F(P l (X),Q i (X,Y))/X Li , 
G l (0,0) = and G iY (0,0) ^ 0. 

By the formal Implicit Function Theorem, the latter conditions ensure that there exists 
a unique power series 5 such that Gi(X, S(X)) = and 5(0) = 0. The corresponding 
parametrization of F is therefore Rj( T) = (Pj(T),Qi(T , S(T) )). The power series 5 can 



be computed using "fast" techniques (|Kung and Traubl (|1978l )'). It may also happen that 



Y divides Gi, in which case the expansion is finite. Therefore, we can consider that such 
a triplet represents a Puiseux series or a rational Puiseux expansion. 

3.2. Classical Puiseux series 



We first give a variant of the classical Newton-Puiseux algorithm (see Walker ( 1978h ) 



to compute the sin gular par t s of P uiseux series. Our floating point algorithm is based 
on this pattern (see Poteaux ( 20071 )). This version returns exactly one representative for 



each cycle above 0. The algorithm is presented in a recursive setting and the output is 
given in a parametric form by mean of triplets as above. The first (non recursive) call to 
the function must be treated differently since £Af(H) must be used instead of QN(H) , for 
reasons explained in Section 4.2. We assume that a mechanism is available to distinguish 
the initial call from recursive calls (adding a Boolean argument would work, for instance). 

CNPuiseux(#) 
Input : 

H : A squarefree polynomial of degree d>\ in L[X, Y], such that H(0,Y) 7^ 0. 
Output : 
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A set of triplets {[Gi, Pi, Qi]}i, which form a set of representatives for : 

- all cycles of H above defined at X = for the initial call, 

- cycles of H that vanish at X = for recursive calls. 

Begin 

If in a recursive call then 

TV <- GN(H) 

If 1(H) = 1 then Return {[H,X,Y]} End 
else 

TV <- £N(H) 
End 

For each side A of TV do 
Compute q, m, I and <pA 
For each (distinct) root £ of 0a do 

a <- f Vfl 

H'(X,Y) <- iJpf«,X m (a + y))/X' 

For each [G, P, Q] in CNPuiseux (i?) do 

K^KU {[G,P«,P m (a + Q)]} 
End 
End 
End 

Return TZ 
End. 

In particular, when applied to the monic polynomial F, CNPuiseux returns a set of 
representatives for all the cycles of F above 0. 

Example 22. Set F(X, Y) = Y 3 -X 5 e Q[X, Y]. Applying CNPuiseux yields one triplet 
with Pi (X) =X 3 andQi(X,F) = X°(0+X 3 (0+X 2 (1+Y))) =X 5 +X 5 Y. The first null 
coefficient comes from the exceptional polygon [(0, 0), (0, 3)]. The second one correspond 
to the fictitious edge of QN(F) introduced at the first recursive call. This may seem 
inefficient, but these tricks have no impact on the complexity and clarifies arguments in 
Section 4. In practice, one may still use a classical Newton polygon if necessary. 

Example 23. Consider F(X,Y) = (Y — 1 — 2X - X 2 )(Y - 1 - 2X - X 7 ) e Q[X,Y]. 
We obtain two triplets with : 

(Pi, Qi) - (X, X°(l + X(2 + X(l + Y))) = (X, 1 + 2X + X 2 + X 2 Y) 
(P 2 , Q 2 ) = (X, X°(l + X(2 + X(0 + Y))) = (X,1 + 2X + X 2 Y) 

Note that the generic Newton polygon allows to obtain immediately the regularity index 
of the series X + X 7 in F, which is 2. The classical polygon does not provide directly 
this information. 

Example 24. Let F be the product of the minimal polynomials over Q(X) of the series 
Jf 5 /6 + X and X 5 ? 6 + X 11 / 12 . We obtain two triplets with : 

(Pi, Qi) = (X 6 , X°(0 + X 5 (l + X(l + Y)))) = (X\ X 5 +X 6 + X 6 Y) 
(P 2 , Q 2 ) = (X 12 , X"(0 + X w (l + X(l + Y)))) = (X 12 , X 10 + X 11 + X^Y) 
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The regularity indices in F are indeed 6 and 11. 



The data q, m, I come directly from the polygons. They need to be computed exactly, 
since for instance q will contribute to the ramification index, which has to be obtained 
exactly. It is obvious that if £ is replaced by a numerical approximation, the change of 
variable in CNPuiseux will almost always produce a polynomial H' with trivial Newton 
polygon, namely reduced to the unique point (0,0). It will not be easy to recover the 
correct polygon, since we will have to decide which coefficients are approximations of 
and should be ignored. Moreover : 

Proposition 25. Let H be a polynomial satisfying the input hypotheses of CNPuiseux. 

• The integer 2(H) is the number of Puiseux series of H above vanishing at X = 0. 

• The integer I(H') is equal to the multiplicity of £ in 4>a- 



Proof. See (iDuvall (|l989h ). □ 



The second assertion of Proposition 25 shows that </>a is not squarefree in general. In 
the presence of approximations, determining the distinct roots of 4>& and their multiplic- 
ity may be difficult. 

However, if we assume that all Newton polygons are obtained by some other means 
(which implies that multiplicities are also known) and provided as an input, then we 
can : 

(1) Extract the approximate coefficients of H which are meaningful to compute QM(H). 
The coefficients below QAf(H) should be equal to : discard them. 

(2) Deduce an approximate 4>A- 

(3) Find clusters of approximate roots of 0a with the expected multiplicities. 

(4) For each cluster, deduce an approximate value of £, apply the numerical change of 
variable and proceed with the recursive call. 

With this technique, we obtain approximate Puiseux series with correct ramification 
indices, that is correct combinatorial data. Our experiments show that approximations for 
the series coefficients are reasonably accurate, yielding accurate evaluations of algebraic 
functions in the neighborhood of critical points. To compute the exact data that are 
needed, we use computations modulo a well chosen prime number p (see Section 4). A 
method to establish the correspondence between data modulo p and numerical data is 



sketched in iPoteauxl (|2007t ). 



Remark 26. For any irreducible F(X,Y) E L[[X]]\Y], at each step of CNPuiseux, the 
polygon Af has exactly one edge and the characteristic polynomial has a unique root. 
Moreover, the sequence of generic Newton polygons encountered depends only on the 
characteristic terms of the Puiseux series (see Section 2.2), and not on the other terms. 
In this sense, these polygons are truly "generic" , since all F with the same characteristic 
yield the same sequence of generic polygons. 
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3.3. Rational Puiseux expansions 

We now describe an algorithm, due to Duval, to compute singular parts of rational 
Puiseux expansions above 0, that we shall use over finite fields. We need two auxiliary 
algorithms, for which we only provide specifications : 

Factor (L,(p) 

Input : 

L : a field 

<p ■ a univariate polynomial in L[T]. 
Output : 

A set of pairs {((pi, fo)}i so that (pi is irreducible in L[T] and (f> = Y\ i 

Bzout (q,m) 
Input : 

q, m : two positive integers 
Output : 

A pair of integers (u, v) so that uq — mv = 1. If q = 1, enforce v — 0, and u = 1. 

Algorithm RNPuiseux (L,H) 

Input : 

L : A field. 

H : A squarefree polynomial of degree d > 1 in L[X 7 Y], such that i?(0, Y) ^ 0. 
Output : 

A set of triplets {[Gj, Pi,Qi]}i, which form a set of representatives for : 

- L-rational Puiseux expansions of H above defined at T = for the 
initial call, 

- L-rational Puiseux expansions of H centered at (0, 0) for recursive calls. 

Begin 

If in a recursive call then 

N <- GAf(H) 

If 1(H) = 1 then Return {[H,X,Y]} End 
else 

N <- £N(H) 
End 

For each side A of Af do 
Compute q, m, I and 0a 
(u, v) <— Bzout(q, m) 
For each (/, k) in Factor (L, 0a) do 
£ <— Any root of / 
H'(X,Y) <- H(£ v X*,X m (C + Y))/X l 
For each [G,P,Q] in RNPuiseux (.L(£), H) do 

K^K U {[G,^P«,P m (C + Q)]} 
End 
End 
End 
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Return 1Z 
End. 

Applied to the monic F, we obtain all rational Puiseux expansions above 0. 
The reader may notice that rational and classical algorithms are quite similar. They 
differ only in the w ay roots of polynomials are handled and in the definit i on of H'. 
In (Duval ( 19891 )1. it is suggested that the D5 system ( Delia Dora et al. ( 19851 )1 should 



be used to avoid factorization. In our case, though, since efficient algorithms are known 
for factoring over finite fields, factorization is not a significant issue (see Section 5). 

Example 27. Let F(X,Y) = {Y 2 - 2X 3 ){Y 2 - 2X 2 ){Y 3 - 2X) G Q[X,Y}. Applying 
RNPuiseux over Q yields three expansions : 

(Pi, Qi) = (2X 2 , X°(0 + 2X 2 (0 + X(2 + Y)))) = {2X 2 ,4X 3 + 2X 3 Y) 
(-P2, Qa) = (4A 3 , X°(0 + X{2 + Y))) = {AX 3 , 2X + 2XY) 
(P3, Qa) = (A, A°(0 + A(V2 + Y))) = {X, V2X + XY) 

The first two expansions have residue field Q and ramification index 2 and 3. The third 
one corresponds to a place with residue field isomorphic to Q(Vz). Applying RNPuiseux 
over Q(V2) will result in one more expansion : 

{Pa, Qa) = (A, A°(0 + X{-^2 + Y))) = (A, - V2A + XY). 

Remark 28. It is worth noting that, unlike classical Puiseux series, rational Puiseux 
expansions are not canonically defined. Replacing T by 7T in Ri{T) = {Xi{T),Yi{T)) 
with 7 chosen in the coefficient field of Ri yields another rational Puiseux expansion 
corresponding to the same place. The choice of 7$ can have dramatic consequences on 
the size of the expansion coefficients and on the performances of the algorithm. See also 
Section 3.6. 

3.4- Polygon trees 

To a function call RNPuiseux(i, F) (see Section 3), we associate a labeled rooted tree. 
By definition, the depth of a vertex v is the number of edges on the path from the root to 
v. In particular, the root vertex has depth 0. The tree vertices of even depth are labeled 
with polygons, while vertices of odd depth are labeled with integer partitions. Similarly, 
tree edges are labeled alternatively with edges of polygons and integer pairs {k, /) where 
k is the multiplicity of a root £ and / = : L}. A tree edge corresponds to the choice 

of a polygon edge or to the choice of a characteristic polynomial root. More precisely, 
the tree is constructed recursively from the root vertex as follow (see Figure 3.4). Even 
depth vertices correspond to function calls. 

• A vertex v of even depth I is labeled with the polygon AT, that is £AT{H) for the root 
vertex (2 = 0), and QN{H) for recursive calls (/ > 0). 

• To each A of Af corresponds an edge from v to a depth I + 1 vertex. Label the edge 
with A (represented by its endpoint). 

• A child (depth I + 1 vertex) is labeled with the corresponding integer partition [4>a] 
(see the end of Section 1 for this notation) . 
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• To each choice of root £ of 0a made by the algorithm corresponds an edge from a 
depth I + 1 vertex to a depth I + 2 vertex. The edge is labeled with the pair (k, /), 
where k is the multiplicity of £ of and / = [£(£) : L]. 

• Then, we proceed recursively : A depth I + 2 vertex is the root vertex of the tree 
corresponding to the function call RNPuiseux(L(£), H') where H 1 is the polynomial 
H' obtained for a choice of edge A and a choice of root £. 



- ((0,0), (7,0)) 

A = ((0,0), (7,0)) 

(7) 

(7,1) 

P = ((0,5), (4,1), (7,0)) 
A = ((4,1), (7,0))/ \a = ((0,5), (4,1)) 

(1) (2, I 2 ) 

(i,i) (2,i) / \a,2) 

V„ V = ((0,1), (2,0)) V„ 
A = ((0,1), (2,0)) 



(1) 



(1,1) 



' = ((0,0), (7,0)) 

A = ((0,0), (7,0)) 

(7) 

(7,1) 

= ((0,5), (4,1), (7,0)) 



A = ((4,1), (7,0))/ 

(1) 

> 1 



A = ((0,5), (4,1)) 



(2, 1 2 ) 



■ ((0,1), (2,0)) V h V h 
A = ((0,1), (2,0)) 



(i) 



Fig. 3. Polygon trees 1ZT(Q, F) and T(F) for Example 27. 



The leaves are even depth vertices labeled with polygons that have only one side 
T^h = [(0, 1), (1, 0)]. Note that the roots £ are not part of the tree. Since the square- 
free factorization is a sub-product of the factorization over L, the labeled tree can be 
obtained at no significant cost. If I is the depth of the function call tree generated by 
RNPuiseux(L,F), then the labeled tree constructed has depth 21. 

For a function call CNPuiseux(F), we define a similar tree, but in this case, an edge 
from a partition to a polygon is only labeled with a multiplicity k because the ground 
field is L and all field extension have degree 1. 

Definition 29. We denote by 1ZT(L, F) (resp. T{F)) a tree associated to the function 
call RNPuiseux(L,i ;l ) (resp. CNPuiseux(F)). In both cases, the tree is called the polygon 
tree associated to the function call. 

It turns out that T{F) is precisely the symbolic information required to achieve our 
goal. 

Proposition 30. The treeT(F) can easily be obtained from 1ZT \L , F) as follow : dupli- 
cate f times each edge labeled (k, f) (together with the sub-tree rooted at this edge) and 
replace the tag (k,f) by the tag k. 
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Proof. From the introduction of Section 3.5, T(F) = 1ZT(L, F). Then, the construction 
is trivial. □ 



This process is illustrated in Figure 3.4. 
3.5. From classical Puiseux series to rational Puiseux expansions 

Following ( Duval ( 19891 ) ). we remark that Newton polygons and root multiplicities 



obtained along the computation are the same with the two algorithms above. This remark 
easily extends to generic Newton polygons and associated characteristic polynomials. 
However, in general, the value of the nonzero roots of characteristic polynomials differ. 

Studying the relations between coefficients of rational Puiseux expansions and classical 
coefficients has a number of benefits : it provides a better understanding of the rational 
algorithm, insight about the coefficient growth (see Section 3.6) and a reduction criterion 
for rational Puiseux expansions (see Corollary 42). Although this criterion is not strictly 
necessary for our algorithm, we have decided to include this material in this paper. 

Let (oci, mi, qi) ■ ■ ■ (cth, rrih, qu) be the sequence of triplets encountered in the compu- 
tation of one classical Puiseux series using the CNPuiseux algorithm. Namely, on is a 
(ft-th root of the i-the characteristic polynomial and —rrii/qi is the slope of an edge of 
the corresponding generic Newton polygon. The output of the algorithm is : 

P{X) = x^-i- = x e 

Q(X, Y) = X m ^-i»(oti + X m2<?3 -«'* (a 2 + • • • + X m »(a h + Y) ■ ■ ■ )) 
so that an element of the corresponding cycle can be written : 

S(X) = X^(ai +X^k{a 2 -\ |-X (a h -\ )•••)) ( 8 ) 

On the other hand, let (£i, mi, qi) . . . (£ r , m/,, qn) be the sequence of triplets encoun- 
tered in the computation of one rational Puiseux expansion using RNPuiseux. Now, £j is 
a root of the z-the characteristic polynomial and —mi/qi is the slope of an edge of the 
corresponding generic Newton Polygon. We denote by (ui,Vi), 1 < i < h the pairs of 
integers returned by the Bzout algorithm. 

Since we have used generic and exceptional Newton polygons, some of the £j and on 
may be null. If = ccj = 0, we have Vi = because is associated with an edge of slope 
-1 or and therefore, % = 1 (see procedure Bzout). In the sequel, we define 0° = 1 so 
that = = 1 and all expressions involved make sense and are correct. 

Proposition 31. There exists a classical Puiseux series as above and a set of integers 
{e-ij}\<j<i<h such that : 



i> = at n 



i-1 

Vi e 



'7 <^17 

a? \ 



Proof. We set Xq = X and Yq = Y, and consider transformations performed by the 
algorithm : 

r,_i = xr(Z?+Y) 
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Wc define (any choice of e-th root is acceptable) : 



so that, we can write : 

X t = mX^-n . 

The truncated series computed by the algorithm can expressed as follow : 

q(x, o) = xr (er + xp (6 + x^ + ■■■ K-x 1 (C-i + x h h e h h )•••)))• 

Using the above expression for Xi and identifying coefficients with those of expression 
(8), there exists a classical Puiseux series so that : 

a i a i _ 1 = e, i _ l VMi l<i<h 

where we have chosen Qq = £o = 1- It is convenient to introduce ft — o^O!^-^"^ 1 ■ Hence, 
we have : 

0i = svr 

Raising the second equality to the power and applying the relation — rriiVi = 1, 
we obtain : 

f. — f) qi ir mi 

Raising the first equality to the power Ui and the second one to the power Vi, Bzout 
relation gives : 

Vi = ft ><-!• 

The recurrence given by the last equality easily yields : 

Q — Vi n — Vi-lUi n~ Vi-2Ui — lUi n — VlU2U3-- Ui 

Vi = ft ft-1 ft-2 , g x 
Si = ft lft-1 ft-2 ■ • - y i J 

Finally, the proposition is proved by induction on i, together with the following assertion : 
There exists a set of integers {fij}i<j<i<h such that : 

i-l 

ft > 1[ a v / fi > (10) 
The case i = 1 is trivial. Assume that i > 1. The induction hypothesis about gives : 

Setting = ei-ijUi-i for 1 < j < i — 2 and fi,j-i = r^i-i wc obtain (10). The 
expression for £j in the proposition then follows directly from the formula for £j in (9) . □ 

Remark 32. Assuming that the Vi are chosen in N, it is easily seen that the and fa 
are also in N. 
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Remark 33. Using relation (9) and the definition of 9i, it is easy to express recursively 
the £j in terms of the on , but there is no simple formula. On the other hand, the on can 
be easily expressed as follow : For 1 < i < j < h, define s ?i = „ m * . Then : 



VjSji 

I 

3 = 1 



To conclude this part, we rewrite the coefficients of parametrizations returned by 
RNPuiseux in terms of the In order to simplify expressions, we introduce the following 
notation for < i < h — 1 : 

( _ t"i+lt1'i+29»+l/ r fi+39i+l9i+2 / r v h q i + 1 ...q h _ 1 
— Si+1 Ci+2 Si+3 ""S/i 

We also define £(m = 1. The output is now given by : 

Q(X, Y) = ^,'17 + ZffixZk (Q 2 +■■■ + + 50 • • • ))■ 

We deduce the parametrization : 

X(T) = e (0) T e 

Y(T) — £" 1 ,£™jT" il92 "' 9h + 

^2 2 ^(1 { ^(2) 2 T 1 " 1,1 92 '" qh+rrl2q3 "' qh + (11) 



5.6. Rational Puiseux expansion coefficient bit-size 

In (|Walshl (|l999h ). it is mentioned without evidence that, unlike classical Puiseux 



series, rational Puiseux expansions can have coefficients with exponential bit-size. Since 
we shall compute rational Puiseux expansions only over finite fields, the coefficient bit- 
size is essentially bounded by the logarithm of the residue field cardinal. Therefore, this 
bad behavior is not a significant drawback for us. 

Example 34. We illustrate the coefficient growth of rational Puiseux expansions and 
the impact of the choice of u and v in the procedure Bzout with the following example. 
The growth is not exponential, but still, leads to dramatically large coefficients. Let a 
and h be positiv e integers and con s ider t he following parametrization, introduced in a 
different context (He nry and Merle (| 19871^ : 



X(T) = T 2 \Y(T) = ^ aT 32 h d-i/2 fc ) (12) 

k=l 

We define d = 2 h and F d (X,Y) = Resu\tenit T (X - X(T),Y -Y(T)), so that degy- (F d ) = 
d and deg x (Fd) = 3(d— 1). There is a unique place above for Fd and therefore, a system 
of rational Puiseux expansion contains a unique parametrization. We use the notations 
of the previous subsection. Applying RNPuiseux to Fd, we see that h is the number of 
recursive calls and mi = 3, q t — 2 for 1 < i < h. 
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• First strategy. We choose Vi nonnegative and minimal, that is Ui — 2 and Vi — 1. 

We have 9\ = a and, for i > 1, 9i = ata~} 1 ^f_ 1 = £ 2 _ 1 so that relation (9) yields : 

6=^i(^t 3 ---er 2 a 21 " 2 ) 3 . 

Therefore, £j is an integer power of a and crude estimates show that the exponent 
is greater than 4 1 for i > 3. Substituting in £( ), we deduce that £( ) is a power of 
a with exponent larger than 4 h 2 h ~ 1 > d 3 /2 for h > 3. Hence, the bit-size of £( ) is 
(much) larger than d 3 /21og(a) for h > 3. The behavior is clearly worse for the other 
coefficients of (11). For F 16 , RNPuiseux returns : 

X(T)=a 3a72 T 16 

Y(T) = a i609 T 2i + a 6913 T 36 fl 80652n42 _|_ ^641/^45 

Considering (12), this result is far from optimal ! 

• Second strategy. Another reasonable choice is Ui = v% = — 1. Hopefully, this will 
result in smaller coefficients since the exponents in (9) will be considerably smaller 
and have alternate signs. For F\q, we obtain : 

X(T)=a 528 T 16 

Y(T) = fl 793 T 24 + a 1189 T 36 + a 1387 T 42 + a 1486 T 45 

Again, we fall short from a satisfactory result. 

• Maple 10 strategy. The algcurves [puiseux] command returns an even worse re- 
sult : 

A(T)=a 24672 T 16 

Y(T) — a 370a9rp2i _|_ a 55513 ji36 _j_ fl 64765y42 a 6939l2n45 



As remarked in iDuvall (|l989l) . optimal output is not always reachable by the current 



algorithm, no matter how u and v are chosen, even for simple cases. Transformations 
different than that of RNPuiseux may be necessary. In this example, reduction of powers 
of a in the course of the computation by substitutions of the form T — U / a s may result 
in smaller coefficients, but it is not clear how efficient this workaround may be in general. 

If L is a number field, Walsh has shown that a coefficient 7 with bit-size in 0{d 11 n 3 ) 
can be construc ted so that X (T) — ^T e and all other expansion coefficients have poly- 
nomial bit-size dWalshl (|l999ft y Walsh's construction requires that a Puiseux series is 



already known and does not allow to optimize the size in the course of the algorithm. It 
is not clear how far this construction is from an optimal size. 

In summary, an algorithm that computes rational Puiseux expansions with provably 
small coefficient size is still unknown. 



4. Good reduction 

We consider a polynomial F with coefficients in an algebraic number field K and 
discuss how to choose a prime number p so that the computation of rational Puiseux 
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expansions modulo p provides enough information to guide floating point computations 
of Puiseux series, namely T(F), 

We denote by o the ring of algebraic integers of K. If p is a prime ideal of o, v p the 
corresponding valuation of K. Finally, we define : 

o p = {a e K | v p {a) > 0}. 
Let L be the finite extension generated over K by the Puiseux series coefficients of 
F, Note that by Proposition 31, L also contains the coefficients of rational Puiseux 
expansions computed by RNPuiseux. If D stands for the ring of algebraic integers of L 
and *P a prime ideal of D, we introduce : 

D<p = {a e L\v<$(a) > 0}. 

In the sequel, will always denote a prime ideal of D dividing p. 

The reduction modulo of a <E D<p is represented by a. We extend this notation 
to polynomials and fractional power series with coefficients in Dtp. If a G o p , since 
divides p, reduction modulo *p and p coincide and we shall use the same notation a. 

4-1. Modular reduction of Puiseux series 

Our reduction strategy is based on the following definition : 

Definition 35. Let p be a prime number and p a prime ideal of o dividing p. If the 
conditions below are verified : 

• Feo p [X,Y], 

• p > d, 

• u p (tc(A F )) = 0. 

we say that F has good reduction at p. 

Note that if F has good reduction at p, since divides p, wqj(tc(A^)) = and F has 
also good reduction at *p. We shall use this fact freely in the sequel. 

A fundamental result for the reduction strategy is the following consequence of a 
theorem by Dwork and Robba : 

Theorem 36. If F has good reduction at p ; then the Puiseux series coefficients of F 
above are in Osp. 



Proof. Let S(X) = J2Zo a * xi/e be an y of the S *j( x )- We write C p for the field of 
p-adic numbers, that is the completion of the algebraic closure of Q p (completion of Q 
for p-adic absolute value). Consider the field L as a subfield of C p by means of its ^-adic 
completion. We denote by |.| p the absolute value of C p , so that : 

Dqj = {ae L\ \a\ p < 1}. 

Since the coefficients of Ap are in Dqj and \tc(Ap )\ p = 1, Ap has no nonzero root i n 
D(0, 1") = {aeC p | H P < 1}. Therefore, Theorem 2.1 of (|Dwork and Robbal (|l979t )) 



asserts that the p-adic convergence radius of S(X) is at least 1. Moreover, S(X) is 
bounded by 1 on D(0, 1~), otherwise, S(x) cannot satisfy a monic equation S(x) d + 
Aa-^Six)^ 1 + ■■■ + A (x) = for x e D(0, 1~). 
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A classical isometry of p-adic analysis (see Robert ( 2000() . Section 4.6) 



sup \S(x)\ p =sup|a i | J) 

igD(0,l-) i>0 



yields \cti\ p < 1 so that on £ Dtp. □ 



It is worth insisting on the fact that this result holds for any dividing p. 

Example 37. Consider the case F(X, Y) = Y 2 — X 3 (p + X) with p > 2. Puiseux series 
above are : 

s^x) = (-i)^ 3/2 (i + -) 1/2 = (-iy^\i + f f- 2 + ■ ■ • ). 

p 2p %p l 

They are obviously not reducible modulo p, but the discriminant criteria of the theorem 
detects this deficiency. It is interesting to note, however, that a system of rational Puiseux 
expansions is given by {X = pT 2 , Y = p 2 T 3 + ^p 2 T 5 + ...}. This parametrization is 
reducible modulo p, but the reduction {X = 0, Y — 0} is trivial and hardly useful. On 
the other hand, {X — T 2 /p, Y — T 3 /p + ^T 5 /p 3 + . . . } is also a (non reducible) system 
of rational Puiseux expansions. 



Remark 38. This reduction criterion could also be deduced from (Dwork (1984)). Under 



the hypotheses of Theorem 36, it is shown therein that Puiseux series modulo p can 
be lifted to Puiseux series with coefficients in C p . The latter must coincide with the 
embedding in C p of Puiseux series over L. Therefore, Puiseux series over L can be reduced 
modulo a prime ideal dividing p. 

Corollary 39. If F has good reduction at p, then any monic factor G of F in K[[X]][Y] 
satisfies wq3(tc(Ac)) = 0. 



Proof. Let F = GH . Both G and H are monic, so relation (2) shows that tc(Af) = 
tc(A(3)tc(Aff)tc(Resultant(G, H)). From Theorem 36, the coefficients of G and H are in 
D<p, and so are these numbers. The result is then a trivial consequence of i> p (tc(Ai?)) = 
u«p(tc(A F )) = 0. □ 



Corollary 40. If F has good reduction atp, then characteristic coefficients of all ramified 
Puiseux series of F above have ty-adic valuation zero. In other words, reduction modulo 
*P preserves the characteristic of ramified cycles of F above 0. 



Proof. Follows easily from Corollary 39 applied to the irreducible factors of F in K [[X]] [Y] 
and from Proposition 6. □ 



It is important to note, however, that annihilation modulo *}3 of Puiseux series coef- 
ficients is not totally controlled by our good reduction criterion. If F is irreducible in 
A'[[X]][T], all non-characteristic coefficients may vanish modulo as shown by Propo- 
sition 6 (consider for instance the minimal polynomial over Q(X) of S(X) = pX + X 3 / 2 ). 
If F is not irreducible, our criterion will also detect cancellation of coefficients that "sep- 
arate" cycles. This property is contained in Theorem 43. 
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Theorem 41. Let {Si}i<j< s be a set of representatives for the cycles of F above 0. 
Assume that F has good reduction at p. Then, {Si}i<i< s form a set of representatives 
for the cycles of F above 0. 

Proof. The ~S~i satisfy F(X,~Sl) = and A F = ^ 0. Therefore, the S~, are distinct 
roots of F. By Corollary 40, the ramification index of Sj is equal to the ramification index 
of Si, namely a. Since X^i=i e » = ^> we have obtained a complete set of representatives 
for the cycles of F. □ 

We now show that parametrizations computed by RNPuiseux yield meaningful parametriza- 
tions when reduced modulo ^3 (see Example 37). 

Corollary 42. Denote by R{T) = {-/T e , £)[ =0 (3 t T a *) (with (3 t ^ 0) a parametrization 
given by RNPuiseux. If F has good reduction at p, then the [3i belong to D<p and vy(j) = 
0. 

Proof. We use the notations of Section 3.5. If 07 is a characteristic coefficient, then 
Corollary 40 shows that v<p(aj) = 0. If aij is not a characteristic coefficient, it is the 
root of a characteristic polynomial of a Newton polygon edge with integer slope. Hence, 
qj = 1, Vj = and Uj = 1 (see procedure Bzout). From Proposition 31 we deduce that 
vrp({;i) — qiV<p(ai) for all 1 < i < h. The same argument proves that wp(£(i)) = for 
< i < h. In particular, 7 = £( ) and ^(7) = 0. Finally, (11) shows that v<p(/3i) = 
UiqiV<p(ai). If cti is a characteristic coefficient, the latter valuation is zero, otherwise it is 
equal to qiV^{cti) > since Ui = 1 in this case. □ 

4-2. Modular reduction of polygon trees 

If F £ o p [X, Y] and p > d, algorithms of Section 3 can be applied to the reduction F 
of F modulo p, so that the notations T(F) and !ZT(¥ p t, F) make sense. The computed 
expansions have coefficients in a finite extension of F p . 

The following result is crucial. It allows to compute by means of modular computations 
the symbolic information required : 

Theorem 43. If F has good reduction at p, then T(F) = T(F). 

Note that the correspondence between T{F) and T(F) cannot be stated so simply if 
classical Newton polygons are used instead of generic ones : Non-characteristic coefficients 
of Puiseux series may vanish upon modular reduction, yielding polygon modifications. 

We begin with a number of lemmas : 

Lemma 44. Assume that H satisfies : 

(i) H£D V [X,Y], 

(ii) H has no multiple roots, deg y H = d > 0, H(0, 0) = 0, H (0, Y) ^ 0, 

(iii) the leading coefficient of H as a polynomial in Y is a power of X, 

(iv) the roots of H are in U e >o£)q3 {{X 1 ^)), 

(v) H has no multiple roots. 

(vi) v v {tc(A H )) = 0. 
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Let (m,q,l) be integers associated to an edge A of QAf(H) and let £ be a root of <j)\. 
Then, H'(X,Y) = H(X<*, X m (£ + Y))/X l also satisfies the conditions (i) to (vi). 

Proof. Conditions H'(0, 0) = and H'(0, Y) ^ follow from the properties of CNPuiseux. 
The leading coefficient property is obviously preserved by the transformation. If {Yi(X)}\< 
denotes the roots of H, the roots of H' are {Yi(X)/X m — £}i<i<d. It is clear that they are 
distinct. Since £ is in £)<p, so are the coefficients of H' and of its roots. Finally, since the 
discriminant is, up to a power of X, a product of root differences, its trailing coefficient 
does not change under the transformation. □ 

Lemma 45. Assume that H satisfies the condition of Lemma 44- Then : 

(i) QN{H) = gN{H). ^ _ 

(ii) Let A be a edge of QAf(H). The characteristic polynomials 4>a (resp. 4>a) of A in 
H (resp. H) satisfy [4>a] = [<M] (equality of root multiplicities). 

Proof. Denote {Si}i<i< w the cycles of H that vanish at X = and {-ffi}i<i<t« their 
minimal polynomials over K((X)). The assumption about the roots of H induces that 
Hi belongs to D<p[[A]][F]. We define V — Y\^ =1 Hi, so that V is a monic polynomial 
with coefficient in Dqj. Define U so that H = UV. 

By Proposition 25, 1(H) = 1(V), thus U(0,0) ^ 0. Necessarily, QM(H) = QM(V). 
We first show that GJ\f(H) = GN(V), which is equivalent to u«p(I7(0,0)) = 0. Since the 
leading coefficient of H is a power of X, relation (2) shows that : 

tc(Ajj) = tc(Ay)tc(A c/ )tc(Resultant(C/, V) f . 

But V is monic, so the latter resultant is ± Y\ { U(X, Vi), where Vi runs through the roots 
of V. Since Vi(0) = 0, it is obvious that tc(A#) is the product of a power of 1/(0,0) 
and of an element of D<$. The hypothesis Ufp(tc(Ajy)) = yields vys(U(0, 0)) = 0. Thus, 

gM(H) = gM(y). 

To prove (i), by Lemma 21, it remains to show that QAf(Hi) = QAf(Hi). If l(Hi) = 1, 
then l(Hi) = 1 because Hi is monic. Therefore, both QAf(Hi) and QJ\f(H) are reduced 
to the unique edge [(0, 1), (1,0)]. Assume that l(Hi) > 1. It is easily seen that the 
hypotheses imply i)(p(tc(Ajj i )) = 0. Proposition 6 shows that Si and 5, have the same 
characteristic. In particular, Hi is irreducible in F P [[X]][Y] and QAf(Hi) has a single 
edge, whose characteristic polynomial has a unique root. If the unique edge of QAf(H^) 
has slope —1, so does the unique edge of QJ\f(Hj) since the vanishing modulo *p of 
tc(Hi(X, 0)) leads to the same (fictitious) edge. If the unique edge has a slope greater 
than — 1, tc(H(X, 0)) is a nonnegative power of a characteristic coefficient and does not 
vanish modulo *p. In both cases, QJ\f(Hj) — QJ\f(Hi). 

To address (ii), let A be a common edge of QJ\f(H) and QJ\f(H). If A corresponds to 
unique irreducible Hi and Hi, <f>\ and 4>a have a unique root with the same multiplicity, 
since they have the same degree, and we are done. Assume that A corresponds to a least 
two irreducible polynomials H x and H 2 associated to the roots £i and £ 2 of 4>a- In order 
to demonstrate (ii), we just need to show that if £i ^ £2, then £1 ^ £2- If rn and q are 
relatively prime integers such that —m/q is the slope of A, we set = ^ q (any choice 
of g-root suits). The cycle associated to Hi can be represented by the series oaX m / q + • • • . 
By (1), there exists 5 G Dqj such that tc(A H ) = (ai — a 2 )5. Hence, v<p(ai — a 2 ) = 0, so 
that en 7^ 02, which in turn gives ^1 / ^- □ 
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Lemma 46. Assume that F has good reduction at p. Then : 

(i) £M(F) — £Af(F) . 

(ii) Let A be the unique edge of £Af(F). The characteristic polynomial 0a (resp. <j>&) 
of A in F (resp. F) satisfy [(f) a] — [<M]- 

Proof. Since F is monic, (i) is trivially correct by definition of £J\f(F). As for (ii), remark 
that A = F(0, T) (resp. Ja~ = F(0,T)) and th at the d roots of r(0, T) (resp. F(0,T)) 
are precisely the d values of Sij(O) (resp. 5^(0)). The good reduction hypothesis implies 
by (1) that if Sij(0) ^ S uv (0), then 5^(0) 7^ 5„„(0). Therefore, the multiplicities of the 
roots of (f>A are preserved by reduction modulo *p. A fortiori, the squarefree factorization 
of (f>A is preserved modulo p. □ 

Remark 47. The assertion (i) of Lemma 46 does not hold if the exceptional polygon is 
replaced by the generic polygon, as shown by the example F(X, Y) = (Y + p + X){Y + 
1 + X). The good reduction criterion does not detect the cancellation of F(Q, 0), but does 
detect a change of root multiplicities. This remark justifies the introduction of £N(F). 

Proof, (of Theorem 43) By Lemma 46, the root vertex, depth one vertices and edges 
down to depth 2 vertices of T(F) are correctly labelled. 

Let £ be a root of F(0, Y) and H(X,Y) = F(X, Y + £). Assumptions of Lemma 44 are 
obviously satisfied for H because £ € £)<p and Ah = Af- Denote by T'(H) the sub-tree 
of T(F) corresponding to the recursive function call CNPuiseux(TJ). 

We show that, for all H satisfying hypotheses of 44, T'(H) = T'(H). We proceed by 
induction on the number c of function calls to CNPuiseux necessary to compute T'(H). 

If c = 1, then T(H) = 1 and T'{H) is reduced to a single vertex labelled with QJ\[(H), 
which consists of the unique edge [(0, 1), (1,0)]. Lemma 45 gives T(H) = T(H). 

Suppose now that c > 1. Lemma 45 shows that the root vertex, the depth 1 vertices and 
all edges from the root to depth 2 vertices of T'(H) and T'(H) coincide and are labelled 
identically Let H' denote a polynomial obtained from H in CNPuiseux. The number of 
function calls necessary to compute T'(H') is less than c. Lemma 44 ensures that the 
induction hypotheses can be applied to H' . Hence, T'(H') = T'(H'). By construction of 
polygon trees, T'(H) = T'(H) □ 

4-3. Choosing a good prime 

This part is devoted to the choice of a prime ideal p such that F has good reduction 
at p. 

Assume that K = Q(j) and let M 7 be the minimal polynomial of 7 over Q. The 
elements of K are represented as polynomials in 7 of degree less than w = [K : Q] with 
coefficients in Q. Up to a change of variable in M 7 and the coefficients of F, we suppose 
that 7 belongs to 0, namely M 7 e Z[T]. 

Definition 48. Let P be a multivariate polynomial of K[X]. There exist a unique pair 
(H,c) with H e Z[T,X], c e N, deg T < w and P(X) = H(j,X)/c, where c is minimal. 
The polynomial H is called the numerator of P and is denoted num(P). The integer c is 
called the denominator of H and is written denom(P). We define the size of P as follow : 
ht(P) =max{logc,log||ii|| 00 }. 
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Defining F n — num(F) and b = denom(_F), we have : 

F(X,Y)= F ^> X > Y \ 
b 

We are left with the problem of finding a prime number p and a prime ideal p of o 
dividing p such that : 
(Ci) p > d. 

(C2) p does not divide b. 

(C3) We can determine an explicit representation of a prime ideal p of dividing p, so 

that a morphism — > o/p = F p * can be effectively computed. 
(C 4 ) tc(A F ) ^ modulo p. 

Condition (Ci) and (C2) are easily verified. We deal with condition (C3) in a standard 
fashion. Let M be any irreducible factor of M 7 in ¥ p [T] and M be a lifting of M in Z[T]. 
It is well-known that if p is a prime number not dividing the index e 7 = [0 : 2(7)] 
and if M is any irre ducible factor of M 7 in F p [T], then the ideal p = (p, M(*y)) of is 
prime (ICohenl Jl993h L Hence, elements of can be reduced by means of the morphism 



— > o/p = F p [T]/(M) = F p t where i = degM. Computing e 7 is a non-trivial task, and 
so is the computation of generators of prime ideals dividing p when p divides e 7 . If e 7 
is unknown, it is sufficient to choose p so that it does not divide Am , since e 7 divides 
Am t . _ 

In practice, M is chosen amongst the factors of Af 7 of smallest degree. Moreover, it is 
worth trying a few primes p in order to reduce t, the case t = 1 being the most favorable. 

As for {C4,), deterministic and randomized strategies are studied in the next subsec- 
tions. In order simplify the analysis, we replace condition (C4) by the following stronger 
condition : 

(C4) Normx/Q(tc(AF)) ^ modulo p. 

If (Ci) to (C4) are verified, then for all prime ideals p dividing p, F has good reduction 
at p. In practice, though, we do not recommend to use (C 4 ). Finally, we introduce the 
notation : 

N F = 6|Norm K/Q (tc(A F ))A A ^|. 
Conditions (Ci) to (C 4 ) are equivalent to : 
(C5) p > d and Np ^ modulo p. 

4-3.1. Deterministic strategy 

We determine a bound B such that, for all prime numbers p > B, condition (C5) is 
satisfied. We prove first two lemmas that will also be useful in the next section. 

Lemma 49. The discriminant A Fn € Z[X, T] of F n with respect to Y satisfies : 

||AfJ|oo < (2d~iy.d d [(w + l)(n+l)] 2d - 2 \\F n \\ 2 ^ 1 . 

Proof. Expanding the discriminant of the Sylvester matrix of F n and F n y , we see that 
there exist d— 1 coefficients {a.i}i<i<d-i of F n and d coefficients {ai+d-i}i<i<d of F n y/d 
such that : 

2ci— 1 

||AfJ|oc < (2d-l)!d d || J] oitx.ioiu 

fc=l 
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The bound follows recursively from : 

IKclloo < (w + l)(n + 1) | |a,i | |oo | |c| |oo 
for all c(X,T) and the inequality ||cii||oo < ||-Fra||oo- a 

We denote by Rf(T) the numerator of tc(A F ). Note that denom(tc(A F )) is a power 
of b dividing b 2d ~ x . 

Lemma 50. 

II^Hoc < HApJIoodlM^loo + i)(— i)(2d-D— +1 = b q (13) 

|Norm x/Q (^( 7 ))| < (w + lf^^M^B- = B 1 (14) 

e 7 < |A Mt | < w w (w + l^-^WM^-i = B 2 . (15) 

Proof. Since the leading coefficient of F n in Y is a constant, evaluation and resultant 
commute and we have : 

6 2d - 1 A F pf) = A F „ (7iX , y) = A F Jj,X). 

Let d (T) £ Z[T] be the coefficient of X 1 in A Fri (T, X) and c t (T) S Z[T] be the numerator 
of the coefficient of X 1 in A Fn ^ X )- It is clear that Ci(j) — C,(7)- Since M 7 is monic, 
Euclidean division yields Qi £ Z[T] such that C; = QiMj + a. Since deg T Cj < (2d — 
l)(w — 1), one can show that : 

IMloo < IIQIUdlM^U + i)(»-D(a«J-i)-»+i. 

The latter inequality gives (13). Now, from Norm K /^(Rp^)) = Resultant (M, Rf), 
Hadamard's inequality and trivial comparison of norms yield (14). The same arguments 
show that (15) holds. □ 

Finally, we have the following result, for which we do not claim optimality : 

Proposition 51. Define B = max{£>, d + l,i?i,£?2} (see (14) and (15)). Then, for all 
p > B condition (C5) is verified. Moreover, B is effectively computable and there exists 
a prime p> B with size ht(p) £ 0(wd(w ht(M 7 ) + ht(F) + log (wind))). 

Proof. If p is a prime greater than B, it is obvious that (C5) is verified. For d > 1, we 
have B\ > d Taking logarithms and using Stirling's formula in the definition of B\ and 
B2, it is readily seen that B has the announced asymptotic size. Since there is always a 
prime between B and 2B, the proposition follows. □ 



4-3.2. Probabilistic strategies 

We begin with a "Monte Carlo" -like method, best described by the following algorithm. 
We need two auxiliary procedures : The function call RandomPrime(^4,C) returns a 
random prime number in the real interval [A, C\. We assume that the n umbers r eturn ed 
are uniformly distributed in the set of primes belonging to [A,C] (see ( Shoup ( 20051 )). 
Section 7.5). The function Nextprime gives the smallest prime larger than the argument. 
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MCGoodPrime (F , M 1 , B' , e) 
Input : 

F : A monic squarefree polynomial of degree d > 1 in K[X, Y]. 
My : A monic irreducible polynomial in Z[T]. 

B' : A real number such that prime factors of Np are less or equal to B' . 
e : A real number with < e < 1. 
Output : 

A prime number p satisfying (C5) with probability no less than 1 — e. 

Begin 

If B' < 100 then Return NextPrime (B') End 
K <- 6 In B 1 1 (e In In B') + 2d/ In d 
C <- mux {2d, K (In K) 2 } 
Return RandomPrime(<i + 1, C) 
End. 



Proposition 52. MCGoodPrime (F ,M lt B' ,t) returns a prime p satisfying 

ht(p) € 0(loglogB' + logd + loge" 1 ). 
In particular, if B' — max {6, B\, B2} (see (14) and (15)), then : 



ht(p) E 0(log (dwlogn) + loght(F) + loght(M 7 ) + loge 
Moreover, the probability that p does not satisfy (C5) is less than e. 

Proof. First note that condition p > d is automatically verified. If n is a positive integer, 
u(n) is classically the number of prime numbers dividing n. For a positive real number 
x, we use the notation ir(x) for the number of primes less or equal to x. Estimates of 
(|Bach and Shallitl (| 19961 ). Section 8.8) give : 



< 7r(x) fx > 17), 7r(cc) < fa: > 1), oj(n) < — (n > 100). 

ma; Im In In n 

For B' > 100, the bound for uj yields : 

to(N F ) < w(denom(F)) + uj(Norm K/Q {R F (T))) + w(Disc(Af~)) < 6 /"^, . 

' In In if' 

If B' < 100, the algorithm returns a correct result of size O(l) with probability 1. We 
therefore assume in the sequel that B' > 100. Let C be a number greater than 2c? and 
r be the probability that a prime given by RandomPrime (d + 1,C) divides N. We just 
have to determine a C large enough so that r < e. But there is always a prime number 
between d and 2d if d > 1, thus 7r(C) — 7r(d) > 1. Since RandomPrime has a uniform 
behavior, we search for a C such that : 

r= ^£) < e . (16) 

7r(C) — 7r(dJ 

Since B' > 100 and d > 1, estimates above show that it is sufficient to find C with : 

6 In B' 2d 



tt(C) >K = 



elnlnB' lnd' 
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Setting C = if(lnif) 2 , we find C/lnC = if (In if ) 2 / (In if + 2 In In if). For £' > 100, if 
is greater than 23 and therefore (lnif ) 2 /(hiif + 2 In In if) > 1. Moreover C > 226 > 17, 
hence : 

"™ *^ K > 

and (16) holds. Moreover, the algorithm returns a prime p with ht(p) = max {log C, log 2c?}. 
Since log C = log if + 2 log log if G 0(log log B' + log d + log e _1 ), the asymptotic result 
follows. For the particular value of B', the asymptotic bound can be deduced from Propo- 
sition 51. □ 



Finally, we consider a "Las Vegas" flavored method : 

LVGoodPrime(F,M 7 ) 
Input : 

F : A monic squarefree polynomial of degree d > 1 in K[X, Y]. 
My : A monic irreducible polynomial in Z[T]. 
Output : 

A prime number p satisfying (C5). 

Begin 

Np <- denom(F)|Norm K/Q ( J R F (T))Disc(M 7 )| 

B' <- max{denom(F), |Norm K/Q (i? F (T))|, |Disc(M 7 )|} 

Repeat 

p <- MCGoodPrime(F, M 7 , B' , 1/2) 
until p does not divide Np End 
Return p 
End. 



Proposition 53. LVGoodPrime(T ; M 7/ ) returns a prime p satisfying : 
ht(p) e 0(log (dio logn) + loght(F) + loght(M 7 )). 
and the average number of iterations is less than 2. 

Proof. Obvious from Proposition 52. □ 

The computation of B' and Np may have a significant cost. In our monodromy context, 
though, we need to determine Ap anyway. Moreover, in practice, we do not compute the 
norm of Ap trailing coefficient, but perform reduction modulo p = (p, M) instead. 

4-4- Computation ofT(F) 

Let if be an algebraic number field and let be the ring of algebraic integers of if. 
We assume that F e K[X, Y] and satisfies the notations and assumptions of section 1. 
If p is a prime ideal of such that F has good reduction at p, the computation of T(F) 
is straightforward : 

• determine a finite field ¥ p t isomorphic to o/p and F, image of F under this isomor- 
phism, 

• compute lZT(F p t , F) using RNPuiseux, 

• deduce T(F) using Proposition 30 (at the cost of a tree traversal), 

• by Theorem 43, T{F) = T{F). 
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5. Complexity of RNPuiseux over a finite field 



In this section, L denotes a finite field and F belongs to L[X, Y}. Otherwise, we keep 
the notations and assumptions of Section 1. We denote by p > d the characteristic of L. 
We also define to = [L : ¥ p ]. 

This section is devoted to the proof of the following theorem : 

Theorem 54. Assuming that FFT-based polynomial multiplication over finite fields 
is used, the RNPuiseux algorithm can compute the singular parts of a system rational 
Puiseux expansions above of F in 0~(d 3 n 2 + d 2 nt$ \ogp) field operations in L. 

As usual the notation CT hides lo garithmic fact ors. 

This result improves the bound of ( Duval ( 19891 )). which is in 0{d 6 n 2 ) field operations. 



Our estimates include factorization cost, while Duval relies on the D5 system to avoid 
factorizations. The gain comes from : 

- truncation of powers of X in the course of the algorithm (see Proposition 57), 

- reducing transformations to shifts of univariate polynomials, for which fast methods 
are available (Proposition 61), 

- a bound for Sf (see below and Proposition 63). 

For our application to the monodromy problem we need to consider expansions above 
all conjugacy classes over L of critical points. More precisely, if Ap = FJi ^ 1S a 
factorization of Ap into a product of irreducible polynomials of £[AT], rational Puiseux 
expansions above roots of are conjugated over L. Therefore, it is sufficient to compute 
a system of rational Puiseux expansions above one root Ci of A^ for each i. We obtain 
the following and remarkably similar theorem : 

Theorem 55. Assuming that FFT-based polynomial multiplication over finite fields 
is used, the RNPuiseux algorithm can compute the singular parts of systems of ratio- 
nal Puiseux expansions above all conjugacy classes over L of critical points of F in 
0~{d?n 2 t^\ogp) field operations in L. 

We first introduce notations and make some assumptions : 

• {Ri}i<i<p with Ri(T) = (Xi(T), Yi(T)) stands for singular parts of a system of rational 
Puiseux expansion above 0, 

• {Gi, Pi, Qi) is the output of RNPuiseux corresponding to Ri, 

• ( r i> e i: /<)) 1 < i < P are respectively the regularity index, the ramification index and 
the coefficient field degree over L of Ri. 

• For each rational Puiseux expansion Ri, we can deduce e^/i Puiseux series denoted 
S ijk (X),l< k<e t ,l<j< fi. 

• Yijk is the singular part of Sijk- 

• We define F = U ijk Y -Y ijk . 

• The following quantity will enter our estimates : 



p 



5f = ^2 f iTi 
i=i 

L t denotes an extension of degree t of L. 
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• M(N) will denote the number of field operations in L t needed to compute the product 
of two polynomials in L t [Z] of degree no larger than N. We recall that M(N) G 0(N 2 ) 
for classical arithmetic and M t (N) G 0(AHog AHoglogiV) c <J~(N) if FFT-based 
multiplication is used. 

• It will be convenient to introduce M(N) = M(N)/N, so that M(N) G O(N) or 
M(N) G 0(logiVloglog7V) C (7(1). 

• A field operation in L t can be made using 0(M(t) logt) field operations in L. 

We refer the reader to ( von zur Gathen and GerhardI (|l999h ) for assertions regarding the 



complexity of operations over finite fields. 

Remark 56. It is worth noting that Sf is an upper bound for the number of elements of 
L necessary to represent the Yi. Indeed, each Yi has at most nonzero coefficients and 
each of those may be represented by at most ft elements of L. Moreover, assume that 
a dense representation is used for the truncated power series Yi (for instance, a vector 
of Ti elements of Lf i ) and assume in turn that the coefficients of Yi are represented by 
vectors of /; elements of L. Then, Sf is precisely equal to the size of the output. 

We split the proof into several results. 

5.1. Truncating powers of X 

We prove the following proposition : 

Proposition 57. Systems of rational Puiseux expansions for F and F Sp above have 
the same singular parts (up to trivial changes of the parameter T). Moreover, singular 
parts of rational Puiseux expansions of F can be computed by applying RNPuiseux to F Sp 
and truncating polynomials H modulo X Sp+1 at each stage of the algorithm. 



Proof. To obtain singular parts of rational Puiseux expansions of F, we just need Gi 
modulo X for all i. Setting Nq — in the bound of Lemma 58 below and applying Lemma 
59, we see that computations modulo Sp + I are sufficient. □ 



To simplify notations in the next lemma, we drop indices : R is a rational Puiseux 
expansion and (r, e, /) and (G, P, Q) are the associated quantities. Let m k i + q k j = l^, 
(1 < k < h) the sequence of Newton polygon edges encountered in RNPuiseux while 
computing (G, P, Q), and F = Hi, . . . , Hh = G be the sequence of input polynomials. 

Lemma 58. Let Nq be a positive integer. In order to compute G N ° , it is sufficient, at 
each stage of the algorithm, to compute modulo X N+1 , where N — -^r + X)fc=i — ~ — ■ 



Proof. The algorithm performs a sequence of substitution : 

Hki^X^iX™*^* + Y)) 



H k+1 (X,Y) = 



X l « 



If H k {X,Y) = Y t! n,,X J Y'. we define H kw (X,Y) = E m , i+%? = w "y^Y*, so that 
H k = ^2 w Hkw Monomials of H kw are transformed into monomials of H k +i arranged 
on the horizontal line j = w — l k (see Figure 5.1). Therefore, to determine G N ° = 
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it is sufficient to know H^ , h >' qh . Proceeding recursively, we obtain the given value of 
N and conclude that it is sufficient to compute all modulo X N+1 . □ 




Fig. 4. Geometric interpretation of a substitution 



This lemma may also be useful if one wishes to compute expansions beyond the reg- 
ularity index. It is clear from the proof that we could adopt an adaptative strategy and 
truncate further by taking into account the information gained at the current stage of 
the execution. It is unlikely that this approach will improve the asymptotic complexity. 

Lemma 59. With the notation of Lemma 58 : Y^i 1 , — - — < Sp. 

J L^ik=\ qi...Q k — r 



Proof. Since F and F have the same rational Puiseux expansion singular parts, apply- 
ing RNPuiseux to F yields the same sequence (m.^, q^, lk)i<k<h &s for F and an output 
(G, P, Q), where P{X) = XX e and Q(X, Y) = Q (X)+X r Y are the polynomials associ- 
ated to F. Set c = Z/j + Ih-iQh + ■ • • + I1Q2 ■ ■ ■ Qh- From the definition of the substitutions, 
we get : 

G(X,Y) = P ^ X \ Q } X ^\ 

Since (1,0) belongs to QN{G), we have vx{Gy(X 7 0)) = 0. Hence, Taylor formula gives 
c = vx{X r F Y (P(X), Qq(X)), which is equivalent to : 

e e e 

Assume that the truncated Puiseux series QodX/X) 1 /") is equal to Y vyk ,. Then : 

J FV(X,Q (pr/A) 1 / e )) = F Y {X,Y Vfk ,) = J] - Y ^)- 



(«,3, k)it{i' ,j',k>) 



Now, since vx^Xvyw - Yijk) < we nave : 



(»,J,*0 1=1 
(i,j,k)^(i> ,j> ,k>) 



□ 
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5.2. Cost of substitutions 



Lemma 60. Let N be a positive integer. Let H £ L t [X, Y], £ G L t , and define H'(X, Y) = 
H(£ v X<>,X m ((, u + Y))/X l , where (m,q,l) corresponds to an edge ofGAf(H). Then, we 

— - N 

can compute the coefficients of H' with 0{N M(d)) field operations in L t . 



Proof. Write H(X, Y) = H W (X, Y) with H w (x, y) = E mi+qj = w ^X^Y\ 

H w {CX\eX rn {l + Y))= Yj a ij(£"X q Y (X m (C +Y)Y 

mi-\-qj—w 

=x w ^ ^cHc + yy 

mi-\-qj—k 

= X W R W (Y + C) 



where R W {Z) = a ijC JJ % % i s a univariate polynomial of degree at most d. 

7ni+qj—w 

We first form the necessary powers of f : The exponent is bounded by (N + l)/q but 
l/q < d since slopes of generic Newton polygons are at least -1. Computing all powers 
up to this bound is achieved in 0{N + d) operations in L t . Constructing the polynomials 
R w requires at most Nd multiplications in L t . Finally, since p > d, the shift in R w 
can be reduced essentially to the multiplication of two degree d polynomials, with cost 
0(M(d)) (|Bini and Pan! (|1994[ )). Since we must perform N such shifts, the total cost is 
in Q(M(d)N). □ 



Proposition 61. The substitutions needed to compute the singular parts of a system of 
rational Puiseux expansion of F requires O \^^M{d)M{d)\ogd\ field operations in L. 

Proof. To compute a rational Puiseux expansion Ri, RNPuiseux performs at most 
substitutions over the field Lf i . Proposition 57 ensures that all substitutions can be 
done modulo X 6p+1 . From Lemma 60, substitutions for computing Ri cost at most 
0{M{d)§pri) operations in Lf i . Taking into account the extension Lf i jL and summing 
over i, we obtain a total cost : 




= 0(5 2 F M{d)M{d) logd) . 

□ 
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5.3. Factorizations cost 

Proposition 62. All factorizations of characteristic polynomials required by RNPuiseux 

can be computed with an expected number of field operations in L : 

0(S F log d [M(d 2 ) + t log pM{d) log d] ) 



Proof. Let L t be the extension of L over which a factorization of 4>a must be determined. 
It is easily seen that degree g?a of 4>a is at most d/t. This is obviously true at the first 
stage of the algorithm, where t = 1. Assume that the property is true at a given stage 
of the algorithm and denote by £ a root of multiplicity k of </>a- If 4>A' is a characteristic 
polynomial of the polygon induced by the choice of £ and d^ is its degree, by Proposition 
25 : 

d A ><k and k[L t (£) : L t ] < d A < d/t. 
Hence, c?a' < d/[X t (£) '■ L] and the property is true by induction. Hence, factorization of a 
characteristic polynomial can be achieved with an average number of 0(M(( j) 2 ) log j + 
tot \ogpM (j) log j) operations in L t ( von zur Gathen and Gerhard! (Il999h . Corollary 



14.30). Each operation in L t can be performed with 0{M{t) \ogt) operations in L. It is 
easily verified that 0(M((f ) 2 ) log f M(t) log*) C 0(M(d 2 ) logd) and 0(M(|) log f) C 
0(M(d) log 2 d), no matter which arithmetic is used. Hence, we obtain : 

O (M(d 2 ) log d + t < log P M (d) log 2 d) . 

Finally, we multiply by r^, bound t by /j and summing over i to obtain the result. □ 



5.4- Bounding 5f 
Proposition 63. 



Sf = / ]nfi < «x(A_f) 



Proof. By definition of the regularity index, for each Puiseux series Sijk(X), there exists 
an other Puiseux series Si j k (X) with io G {1 ... p}, jo € {1 . . . /,„} and fco € {1 . . . e^} 
such that arl < ^(Sy fc (X)-S ioiofco pO) < a. Moreover, if i- x (S'y fc (^)-5* ioJofco (X)) ^ 
I 1 , &i is a proper divisor of ej . Thus, denoting q = e^/e^ > 1, there exists m £ N and 
a^OgL such that 1 < m < q and : 

Si j k (X) = S m H (X) + aX Cl +c 'o -| 

Hence, for < I < q — 1, we have : 

^W^ W = ^(X) + C?'ax'-^ + ^ +■■■ 
Therefore, we obtain : 

E^( s m(x)-s^(x)) = q!l + ^>!l 

i=a €l 6l Gl 
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For each Puiseux series Sijk, we have : 

v x {F Y {X,S ijk (X))) = J2 v x {S ijk {X) - S Vfk ,(X)) > r -± 



(i'.j'.k') 



Summing over (i,j,k), relation (1) gives : 



v x (a f )= ]r MW.^*W))>EEES: = E r ^ 

(i,j,k) i=l j=l k=l 1 i=l 



□ 



5. 5. Proof of Theorem 54 and 55 

It is interesting to bound first the number of operations in L in terms of the output 
size, namely S F . 

Theorem 64. The number of operations in L required to compute singular parts of 
rational Puiseux expansions of F above is in : 

a (S F M(d) [S F + M(d) + t a logp]) C a (n M(d)d [nd + M(d) + t logp]) . 

Proof. The first assertion follows from Proposition 61 and Proposition 62 since M(d 2 ) € 
0(M(d) 2 ) and M(d)\ogd E CT(1). The inclusion is a consequence of Proposition 63, 
where vx(A F ) is bounded from above by deg Ax < (2d — l)n. □ 



Theorem 54 is now a trivial consequence of this result. 
As for Theorem 55, we proceed as follow : 

First of all, the cal culation of the discriminant can b e done in 0(nM(nd) \og(nd)) 
field operations in L (jvon zur Gathen and Gerhard! ^M)). Thus, this step is included 
in our complexity bound. 

Moreover, is a polynomial in X with degree 8 at most (2d — l)n. Therefore, it can 
be factorized over L in 0(M(5 2 ) log<5 + t log pM (S) log 6) C (J(n 2 d 2 + ndtpl ogp) field 
operations in L using fast multiplication ( von zur Gathen and Gerhard ( 19991 ) V Hence, 
this step is also included in our complexity bound. 

Then, if Ap — II^Li ^ s the factorization obtained, set ti = deg x (Ai) and let Ci 
be a root of A,. The computation of Fi = F(X + Ci, Y) can be performed at the cost of 
d shifts in L(ci) — L ti of the coefficie nts of F in Y . Hence , the complexity of this step is 
in 0(dM(n)) field operations in L ti (|Bini and Pan! (|l994b ) and 0(dM(n)M{U) log U) C 
CT (dnti) operations in L. Summing over i and bounding ^ ti by (2d — l)n gives again 
a correct estimate. 

Finally, we note that the quantity S Fi associated to Fi is bounded by ki (see Proposition 
63). From Theorem 64, using fast multiplication, the function call RNPuiseux(F(X + 
Ci, Y), L(ci)) requires CT (kid[ki+d+toti logp]) field operations in L ti , and so (J~ (kitid[ki + 
d + totilogp]) C (T (kitind 2 to logp) field operations in L. Summing over i allows to 
conclude since J^i = ^ — (2d — l) ?l - 
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6. Bit-complexity 

Let F be a polynomial of K[X, Y], where K is an algebraic number field represented 
as in Section 4.3. We recall that [K : Q] = w. We study the bit-complexity of the com- 
putation of T(F). We estimate only word operations generated by arithmetic operations 
in various coefficient fields. Assuming some care is taken in the implementation, (for in- 
stance, access to coefficients of polynomial should be achieved in constant time), results 
below should give a realistic upper bound for the behavior of an actual program. 

Our bounds for randomized algorithms do not include the cost of generating prime 
numbers, nor the cost of computing bounds given by our formula. 

We assume that elements of ¥ p are represented by nonnegative integers. In order to 
simplify expressions, we assume that fast arithmetic is used for integer arithmetic as well 
as polynomial arithmetic over finite fields. 

Theorem 65. Given a real number e with < e < 1, there exists a probabilistic Monte 
Carlo algorithm that computes T{F) with an expected number of : 

<T(dVw 2 ht(p) 2 [ht(M 7 ) + ht(F)]) C Cr(d 3 n 2 w 2 log 2 e _1 [h.t(M 7 ) + ht(F)]) 

word operations and probability of error less than e. 

Proof. Call MCGoodPrime(f 1 ,M 7 ,.B / ,e) with B' = max {b, Bi,B 2 } to get a good prime p 
with probability of error less than e and size ht(p) given by Proposition 52. Reduce the 
coefficients of M 7 and F n modulo p, and check that p does not divide b — denom(F). 
Denote F n the reduction of F n . This step requires (T (wht(M 7 ) + wndht(F n ) + ht(6)) logp 
word operations. Then, factorize M 7 over ¥ p at the cost of 0~ (w 2 + wlogp) operations 
in F p and choose an irreducible factor M of minimal degree. Coefficients of F n must 
be reduced modulo M and p. There arc nd coefficients and each division induces O" (w) 
operation in F p . The worst case for RNPuiseux occurs if M happens to have degree 
w. In this case, call RNPuiseux(F p ™,F), which require CF (d 3 n 2 + d 2 nw\ogp) operations 
in Fpu, by Theorem 54, each operation requiring CT (w) operations in F p . In summary, 
CT (d 2 nw[dn + wlogp + ht(M 7 ) + ht(F„)] \ogp) word operations in are needed. Taking 
into account the bound for ht(p) = logp and applying crude estimates yield the result, 
since loght(F) and loght(M 7 ) are logarithmic terms in the size of the input. □ 

Remark 66. It is possible to demonstrate without much difficulty that the same bit- 
complexity result holds for the computation of polygon trees of F above all critical points. 
In turn, this result gives a bit-complexity bound for the computation of the genus of an 
algebraic curve, a problem for which we have not found similar results in the literature. 

7. Conclusion 

This paper summarizes the results we have obtained regarding the symbolic part of 
our program towards a fast and reliable method to compute Puiseux series with floating 
point coefficients. In particular, the criterion ensuring preservation of polygon trees is 
essential. 

Along the path, we have derived improved complexity bounds for the computation of 
Puiseux series over finite fields. Although not optimal, these bounds are quite reasonable, 
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i.e. quadratic in the output size, up to a factor d. Refinements may be obtained by 
optimizing several steps of our complexity analysis. It is not clear that these refinements 
may significantly improve the final asymptotic cost. 

Bit-complexity estimates for the Monte-Carlo version of the first stage of our symbolic- 
numeric method confirm that the coefficient growth of pure symbolic Newton-Puiseux 
algorithm is avoided. Complexity bounds for the Las Vegas and deterministic versions 
can be obtained similarly. An analogous bit-complexity bound can also be deduced for 
the computation of the genus of an algebraic curve defined over an algebraic number 
field. 

The numerical part of the algorithm is still being investigated and several variants are 
being compared. 

The scope of our symbolic-numeric strategy may be broader : For instance, series 
solutions at singular points of certain linear differential equations can be constructed 
using Newton polygons. It would be interesting to investigate the benefits of our approach 
in this context. 
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